Transverse Instability of Avalanches in Granular Flows down Incline 
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Avalanche experiments on an erodible substrate are treated in the framework of "partial flu- 
idization" model of dense granular flows. The model identifies a family of propagating soliton-like 
avalanches with shape and velocity controlled by the inclination angle and the depth of substrate. 
At high inclination angles the solitons display a transverse instability, followed by coarsening and 
fingering similar to recent experimental observation. A primary cause for the transverse instability is 
directly related to the dependence of soliton velocity on the granular mass trapped in the avalanche. 

PACS numbers: 47.10.+g, 68.08.-p, 68.08.Bc 



Granular deposit instabilities are ubiquitous in na- 
ture; they display solid or fluid-like behavior as well 
as catastrophic events such as avalanches, mud flows 
or land slides. A somewhat similar phenomena un- 
fold below sea level. Their occurrence is relevant for a 
broad variety of marine-based technologies, such as off- 
shore oil exploitation or deep-sea telecommunication ca- 
bles, and is a matter of concern for coastal communi- 
ties. The perspective of risk modelling of these unsta- 
ble matter waves is hindered by the lack of conceptual 
clarity since the conditions triggering avalanches and the 
rhcology of the particulate flows are poorly understood. 
While extensive laboratory-scale experiments on dry and 
submerged granular materials flowing on rough inclined 
plane 0, 0, [j, 0, IM EJ have brought new perspectives for 
the elaboration of reliable constitutive relations, many 
open questions still remain such as avalanches propaga- 
tion on erodible substrates. It has been shown exper- 
imentally that families of localized unstable avalanche 
waves can be triggered in the bi-stability domain of phase 
diagram Also, the shape of localized droplet-like 
waves was recently shown to depend strongly on the 
intimate nature of the granular material used All 
these questions are closely related to the compelling need 
for a reliable description of the fluid/solid transition for 
particulate assemblies in the vicinity of flow arrest. Re- 
cent avalanche experiments on erodible layers performed 
both in air and under water Q though strongly differ- 
ing by spatial and time scales involved, display striking 
common features: solitary quasi one-dimensional waves 
transversally unstable at higher inclination angles. The 
instability further develops into a fingering pattern via a 
coarsening scenario. So far, this phenomenology, likely 
to be common to many natural erosion/deposition pro- 
cesses, misses a clear physical explanation. From a theo- 
retical perspective, a model of "partially fluidized" dense 
granular flows was recently developed to couple a phe- 
nomenological description of a solid/fluid transition with 
hydrodynamic transport equations. It reproduces many 
features found experimentally such as metastability of a 



granular deposit, triangular down- hill and balloon- type 
up-hill avalanches and variety of shear flow instabilities 
Q. The model was later calibrated with molecular 
dynamics simulations 0. 

In this Letter the partial fluidization model is applied 
to avalanches on a thin erodible sediment layer. A set of 
equations describing the dynamics of fully eroding waves 
is derived, and a family of soliton solutions propagating 
downhill is obtained. The velocity and shape selection 
of these solitons is investigated as well as the existence 
of a linear transverse instability. The primary cause of 
the instability is identified with the dependence of soli- 
ton velocity on its trapped mass. A numerical study 
is conducted to follow nonlinear evolution of avalanche 
front. All these features are discussed in the context of 
the experimental findings of Malloggi et al.Q. New per- 
spectives for quantitative contact between modelling and 
experiments are then underlined. 

According to the partial fluidization theory , the ra- 
tio of the static part of shear stress to the fluid part of 
the full stress tensor is controlled by an order parameter 
(OP) p, which is scaled in such a way that in granular 
solid p — 1 and in the fully developed flow (granular liq- 
uid) p — ► 0. At the "microscopic level" OP is defined as 
a fraction of the number of persistent particle contacts to 
the total number of contacts. Due to a strong dissipation 
in dense granular flows, p is assumed to obey purely re- 
laxational dynamics controlled by the Ginzburg-Landau 
equation for generic first order phase transition, 
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Here T p ,l p ~ d are the OP characteristic time and length 
scales, d is the grain size. F(p, S) is a free energy density 
which is postulated to have two local minima at p = 1 
(solid phase) and p = (fluid phase) to account for the 
bistability near the solid-fluid transition. The relative 
stability of the two phases is controlled by the param- 
eter S which in turn is determined by the stress ten- 
sor. The simplest assumption consistent with the Mohr- 



Coulomb yield criterion is to take it as a function o: 
4> = max \cr mn /cr nn \, where the maximum is sought ove] 
all possible orthogonal directions m and n. 

For thin layers on inclined plane Eq. can be 

simplified by fixing the structure of OP in z-directior 
(z perpendicular to the bottom, x is directed dowr 
the chute and y in the vorticity direction) p = 1 - 
A(x,y) sin(irz/2h), h is the local layer thickness, A if 
slowly-varying function. This approximation valid fo] 
thin layers when there is no formation of static layer be- 
neath the avalanche. Then one obtains equations gov- 
erning the evolution h and A, coordinates x, y, height Jh 
and time t are normalized by l p , t p correspondingly 
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p, is the shear viscosity, tp is the chute inclination, <f> — 
tan tp. Control parameter 8 includes a correction due tc 
the change in the local slope 8 — 80 + (3h Xl 3 « 1.5 — ; 
depending on the value of tp, see for detail [ZLul- The lasl 
term in Eq. (0 is also due to change of local slope anc 
is obtained from expansion tp = Cp + h x . This term is re- 
sponsible for the saturation of the slope of the avalanche 
front (without it the front can be arbitrary steep) . 

In the coordinate system co-moving with the velocity 
V Eqs. ©,© assume the form: 
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Numerical studies revealed that the one-dimensional Eqs. 
© , iJBJ possess a one-parametric family of localized (soli- 
tons) solutions, see Fig^ 



A(x, t) = A(x -Vt),h(x,t) =h(x- Vt) 
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Here the boundary conditions take a form h — > ho , A — » 
for x —> ±00, where ho is the asymptotic height. The 
one-dimensional steady state soliton solution Q satisfy: 
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The solutions can be parameterized by the "trapped 
mass" m carried by the soliton, i.e. the area above ho, 
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FIG. 1: h (a) and A (b) for various values of m and a. Solid 
line is for m = 147.7, V = 2.72, dashed line is for m = 211, 
V = 3.12, for 8 = I, a = 0.08, 8 = 2; point-dashed line is for 
a = 0.025, 5 = 1.15, m = 62, V = 0.86. Inset: V vs m. 



The velocity V is an increasing function of m, see inset 
Fig. [Ql. The structure of the solutions is sensitive to the 
value of a: for large a the solution has a well-pronounced 
shock-wave shape, Fig. ^ with the height of the crest 
hmax several times larger than the asymptotic depth ho- 
For a — > the solution assumes more rectangular form, 
see Fig. [TJ and h max - h < ft. - 

To understand transverse instability we focus on the 
soliton solution with slowly varying position xo(y,t) 

A(x, t) = A(x - x {t, y)), h(x, t) = h(x - x {t, y)) (11) 

Substituting Eq. (|llf> in Eq. (jSJ and integrating over x, 
one obtains 

d t m = V{m)(h+ - h-(m)) - dd 2 x + C, 2 d 2 y m (12) 
where (1^ = const is defined as 



(Ah 3 d x h) dx, (2 = - (Ah 3 d m h) dx 



Ci = - 



Here h + = h(x — > 00) is the height of the deposit layer 
ahead of the front and h~ = h(x — * —00) is the height be- 
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FIG. 2: \{q) vs q for S = 1.15, a = 0.08 and m = 102. 
Solid line: X(q) obtained by numerical stability analysis of 
one-dimensional solution Eq. II II . Dashed line is solution of 
Eq. 1151 . Inset: optimal wavenumber of q* vs a for 8 — 1.15 

hind the front, see Fig. ^l. While the value of h + is pre- 
scribed by the initial sediment height, the value of h~ be- 
hind the front is determined by the velocity (or mass) of 
the front. For steady-state solution h + — h~ — ho- For 
the slowly-evolving solution the difference between h + 
and h~ can be small, however it is important for the sta- 
bility analysis. These terms are also necessary to describe 
experimentally observed initial acceleration/slowdown of 
the avalanches. Substituting Eqs. I|ll|l into Eq. and 
performing orthogonality conditions one obtains 

d t x = V(m) + d*x (13) 

There are also higher order terms in Eq. Q13[) which we 
neglect for simplicity. To see the onset of the instability 
we keep only the leading terms in Ea. (|12|) . l|13[l . using 
V(m) « V(mo) + V m (m — mo), and rh — m — too <C mo: 

d t m = -T-m - C,\d 2 yX {) + C2^m 

d t x = V m ih + dyX (14) 

where mo = const is the steady-state mass of the soli- 
ton, and r = V(mo)d m h ~. Seeking solution in the form 
m,XQ ~ exp[At + iqy], q is the transverse modulation 
wavenumber, for the most unstable mode we obtain from 
Eq. the growthrate A 

. -g 2 (l + C2) - r + y/(g 2 (l - C 2 ) - r) 2 + W m Cig 2 

2 

(15) 

Expanding Eq. (f 1 3|) for q — > we obtain A w 
\{2V m Ci/r - l)q 2 + 0{q A ). The instability occurs if 
V m (i/T — 1/2 > 0. Substituting t and using V m /h m = 
Vh, we obtain a simple instability criterion: 

2V h d/V > 1 (16) 

Eq. (|16fl gives a value of threshold a since £i ~ a. For 
a < a c no instability occurs, and the modulation wave- 



length diverges for a — > a c . Far away from the threshold 
we neglect r and then obtain for X(q): 

\=\q\V(^-(l + (2)q 2 /2 + 0(q 3 ) (17) 

The optimal wavenumber q* is given 

q*-VC^-a (18) 

Fig. [21 shows A(g) obtained by numerical stability anal- 
ysis of linearized Eqs. |J2J, ijSJ near the one-dimensional 
solution Eq. J7J). For comparison is shown the solu- 
tion to Eq. (|15fl . w ith the parameters extracted from 
the corresponding one-dimensional steady-state problem 
Eqs. One sees that Eq. fTK|) gives correct de- 

scription for small q, however fails to predict X(q) in the 
whole range of q. For this purpose one needs to include 
higher order terms. Thus, Eq. I|15|) gives a correct de- 
scription of the onset of instability and qualitative es- 
timate for the selected wavenumber q* . Inset to Fig. 
El shows the dependence of optimal wavenumber q* vs 
a, obtained by numerical linear stability analysis of the 
soliton solution. It shows almost linear decrease of q* 
with a consistent with Eq. I|18fl . F° r very small a the 
plot indicates that q* — > at a — > a c , consistent with 
Eq. H16(l . From the qualitative point of view, the trans- 
verse instability of planar front is caused by the following 
mechanism: local increase of soliton mass results in the 
increase of its velocity and, consequently, the "bulging" 
of the front. Due to the mass conservation, the bulge 
depletes material in the neighboring areas and further 
decreases their speed. 

To study the evolution of the avalanche front be- 
yond the initial linear instability regime, a fully two- 
dimensional numerical analysis of Eqs. (J2J, © was per- 
formed. Integration was performed in a rectangular do- 
main with periodic boundary conditions in x and y direc- 
tions. The number of mesh points was up to 1200 x 600 
or higher. As an initial condition we used a flat state 
h = ho with a narrow stipe h = ho + 2 deposited along the 
y-direction. To trigger the transverse instability, small 
noise was added to the initial conditions. The initial con- 
ditions rapidly developed into a quasi-one-dimensional 
solution described by Eq. 0. Due to the periodicity in 
the x-direction, the soliton could pass through the inte- 
gration domain several times. It allowed us to perform 
analysis in a relatively small domain in the x-direction. 
The transverse modulation of the soliton leading front 
was observed after about 100 units of time for the pa- 
rameters of Fig. |31 We observe that modulation initially 
grows in amplitude, eventually coarsens and leads to the 
formation of large-scale finger structures. 

At the qualitative level the agreement between theory 
and experimental results of Mallogi et al. A\ is impres- 
sive, (i) Existence of steady-state soliton-like avalanches 
propagating downhill with a shape similar to experiment. 
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is proportional to the sediment height ho which increases 
as the angle tp decreases, it results in the decrease of r p . 
Thus, with the decrease of angle <p the instability should 
disappear, in agreement with experiment where the soli- 
ton is found stable at lower inclination angles. 




FIG. 3: Grey-coded images of h(x,y) (white corresponds to 
larger h) for a) t = 170, b) t = 300 and c) t = 500 units of 
time. Domain size is 600 units in x and 450 units in y direc- 
tion, only part of domain in x direction is shown. Parameters: 
S = 1.16, a = 0.14, p = 2 and initial height ho = 2.285. 

(ii) Generic zero wave number (longwave) transverse in- 
stability compatible with the experimental divergence of 
the selected wavelength close to the instability threshold. 
Far from the threshold, linear growth rate dependence 
with q compatible with measurements, (iii) Coarsening 
in the later development of the instability, (iv) Fin- 
gering instability with localized droplet-like avalanches 
(also similar to those described in 01). The analysis pre- 
dicts that the transverse instability ceases to exist when 
the rescaled transport coefficient a decreases (see Fig. 
|2J). In the present form, the model does not provide 
an explicit relation between a and the chute angle tp 
(since a depends also on r p ). Nevertheless, molecular dy- 
namics studies indicate that the OP diffusion coefficient 
Dp = lp/Tp increases with pressure |9|- Since the pressure 



An important question remains is how to bring to a 
more quantitative level the comparison between theory 
and the experimental measurements. In this perspec- 
tive, a challenging question is to deeply understand the 
qualitative differences between smooth glass bead and 
rough sandy materials as far as the effective flow rules and 
avalanche shapes are concerned. This work calls for more 
systematic measurements centered on the soliton velocity 
dependence with the flowing mass for various materials 
and the possible identification of an instability threshold 
for glass beads. Such results would allow a more precise 
assessment of the model parameters and could lead the 
way to a reliable and predictable modelling of granular 
avalanches. The fingering patterns bear remarkable sim- 
ilarities with those existing in thin films flowing down 
inclined surfaces, both with clear and particle- laden flu- 
ids ^(J- However, the physical mechanisms leading to 
this fingering are likely dissimilar: in fluid films, it is 
driven (and stabilized) by the surface tension, whereas in 
the granular flow case, the surface tension plays no role. 
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